jCO2 analysis
Ross Cunning
October 5, 2016
Code for new model
The rejection flux of carbon from host biomass SU, \(jCw\), is used as input to photosynthesis SU. A multiplier is included indicating the number of moles of CO2 the host can deliver to the photosynthesis SU using one mole of fixed carbon. Mechanistically, the host will use this carbon to generate ATP to energize a variety of CO2-concentrating mechanisms (CCM’s) that enhance CO2 availability to photosynthesis.
run_coral_wC <- function(time, env, pars, NPQ="inhibited", kCO2=5) {
require(dplyr)
# Define parallel complementary synthesizing unit function
synth <- function(x, y, m) 1/((1/m)+(1/x)+(1/y)-(1/(x+y)))
#synth <- function(x, y, m) (m^-2 + x^-2 + y^-2)^(-1/2)
# Set initial values
# ==================
# Host
H <- data_frame(
time=time,
H=pars$initH,
jX=(pars$jXm * env$X[1] / (env$X[1] + pars$KX)),
jCO2a=pars$jCO2a,
jCO2p=pars$jCO2p,
jCO2=jCO2a+jCO2p,
jN=(pars$jNm * env$N[1] / (env$N[1] + pars$KN)),
rhoN=jN,
jCw=10,
jHG=0.25,
jHT=pars$jHT0,
rNH=jHT * pars$nNH * pars$sigmaNH,
rCH=jHT * pars$sigmaCH,
dH.Hdt=0
)
H[2:nrow(H), 2:ncol(H)] <- NA
# Symbiont
S <- data_frame(
time=time,
S=pars$initS,
jL=env$L[1] * pars$astar,
jCP=max(0, synth(jL * pars$nLC, (H$jCw[1] + H$jCO2p[1] + H$jCO2a[1])*H$H/S, pars$jCPm), na.rm=T),
jeL=max(jL - jCP/pars$nLC, 0),
jNPQ=pars$kNPQ/pars$nLC,
jCO2w=(H$jCw + H$jCO2)*H$H/S - jCP,
jSG=pars$jSGm/10,
rhoC=jCP,
jNw=0,
jST=pars$jST0,
rNS=jST * pars$nNS * pars$sigmaNS,
rCS=jST * pars$sigmaCS,
cROS=1,
dS.Sdt=1)
S[2:nrow(S), 2:ncol(S)] <- NA
# Run simulation by updating
# ==========================
for (t in 2:length(time)) {
# Photosynthesis SU
# =================
# Light input flux
S$jL[t] <- (1.256307 + 1.385969 * exp(-6.479055 * (S$S[t-1]/H$H[t-1]))) * env$L[t] * pars$astar
# CO2 input flux
H$jCO2p[t] <- pars$jCO2p # passive CO2 flux from the environment
H$jCO2a[t] <- pars$jCO2a # active delivery of CO2 from the environment (mediated by the host)
H$jCO2[t] <- H$jCO2p[t] + H$jCO2a[t] # Total CO2 input from the environment
S$rCS[t] <- pars$jST0 * pars$sigmaCS # metabolic CO2 recycled from symbiont biomass turnover
H$rCH[t] <- H$jHT[t-1] * pars$sigmaCH # metabolic CO2 recycled from host biomass turnover
H$jCw[t-1] <- H$jCw[t-1] # carbon rejected from host biomass SU, given back as CO2 to photosynthesis
# Production flux (photosynthetic carbon fixation)
S$jCP[t] <- synth(S$jL[t] * pars$nLC, (kCO2*H$jCw[t-1] + H$jCO2[t] + H$rCH[t])*H$H[t-1]/S$S[t-1] + S$rCS[t], pars$jCPm) / S$cROS[t-1]
# Rejection flux: CO2 (wasted to the environment)
S$jCO2w[t] <- max((H$jCO2[t] + H$rCH[t])*H$H[t-1]/S$S[t-1] + S$rCS[t] - S$jCP[t], 0)
# Rejection flux: excess light energy not quenched by carbon fixation
S$jeL[t] <- max(S$jL[t] - S$jCP[t]/pars$nLC, 0)
# Amount of excess light energy quenched by NPQ
#S$jNPQ[t] <- min(S$jeL[t], (S$jCP[t] * pars$kNPQ)/pars$nLC)
if (NPQ=="inhibited") S$jNPQ[t] <- min(S$jeL[t], pars$kNPQ / S$cROS[t-1])
if (NPQ=="inhibited2") S$jNPQ[t] <- (pars$kNPQ^(-3)+S$jeL[t]^(-3))^(-1/3) / S$cROS[t-1]
if (NPQ=="fixed") S$jNPQ[t] <- min(S$jeL[t], pars$kNPQ)
if (NPQ=="fixed2") S$jNPQ[t] <- (pars$kNPQ^(-3)+S$jeL[t]^(-3))^(-1/3)
# Scaled ROS production due to excess excitation energy (=not quenched by carbon fixation AND NPQ)
S$cROS[t] <- 1 + ((S$jeL[t] - S$jNPQ[t]) / pars$kROS)^pars$k
# Symbiont biomass SU
# ===================
# Nitrogen input flux
S$rNS[t] <- pars$jST0 * pars$nNS * pars$sigmaNS # Recylced N from symbiont biomass turnover.
H$rhoN[t-1] <- H$rhoN[t-1] # Nitrogen shared from the host (defined below, so previous time step used)
# Carbon input flux
S$jCP[t] <- S$jCP[t] # Production of fixed carbon from photosynthesis SU
# Production flux (symbiont biomass formation)
S$jSG[t] <- synth(S$jCP[t], (H$rhoN[t-1]*H$H[t-1]/S$S[t-1] + S$rNS[t])/pars$nNS, pars$jSGm)
# Rejection flux: carbon (surplus carbon shared with the host)
S$rhoC[t] <- max(S$jCP[t] - S$jSG[t], 0)
# Rejection flux: nitrogen (surplus nitrogen wasted to the environment)
S$jNw[t] <- max(H$rhoN[t-1]*H$H[t-1]/S$S[t-1] + S$rNS[t] - pars$nNS * S$jSG[t], 0)
# Symbiont biomass loss (turnover)
S$jST[t] <- pars$jST0 * (1 + pars$b * (S$cROS[t] - 1))
# Host biomass SU
# ===============
# Food input flux (prey=both carbon and nitrogen)
H$jX[t] <- (pars$jXm * env$X[t] / (env$X[t] + pars$KX)) # Prey uptake from the environment
# Nitrogen input flux
H$jN[t] <- (pars$jNm * env$N[t] / (env$N[t] + pars$KN)) # N uptake from the environment
H$rNH[t] <- H$jHT[t-1] * pars$nNH * pars$sigmaNH # Recycled N from host biomass turnover
# Carbon input flux
S$rhoC[t] <- S$rhoC[t] # Carbon shared by the symbiont (defined above)
# Add dissolved organic carbon uptake
# Production flux (host biomass formation)
H$jHG[t] <- synth(S$rhoC[t]*S$S[t-1]/H$H[t-1] + H$jX[t], (H$jN[t] + pars$nNX*H$jX[t] + H$rNH[t]) / pars$nNH, pars$jHGm)
# Rejection flux: nitrogen (surplus nitrogen shared with the symbiont)
H$rhoN[t] <- max(H$jN[t] + pars$nNX * H$jX[t] + H$rNH[t] - pars$nNH * H$jHG[t], 0)
# Rejection flux: carbon -- given back to symbiont as CO2 input to photosynthesis
H$jCw[t] <- max(H$jX[t] + S$rhoC[t]*S$S[t-1]/H$H[t-1] - H$jHG[t], 0)
# Host biomass loss
H$jHT[t] <- pars$jHT0
# State equations
# ===============
# Specific growth rates (Cmol/Cmol/d)
S$dS.Sdt[t] <- S$jSG[t] - S$jST[t]
H$dH.Hdt[t] <- H$jHG[t] - H$jHT[t]
# Biomass (Cmol)
S$S[t] <- S$S[t-1] + S$dS.Sdt[t] * S$S[t-1] * (time[2] - time[1])
H$H[t] <- H$H[t-1] + H$dH.Hdt[t] * H$H[t-1] * (time[2] - time[1])
}
# Return results
# ==============
return(list(time=time, env=env, pars=pars, H=H, S=S))
}Comparison of new and old model behavior
Varying nitrogen
Old model
# Set up simulation
time <- seq(1,1230,1)
defpars <- def_pars()
defpars <- replace(defpars, "jHGm", 1)
env <- init_env(time=time, L=c(15,15,2), N=c(1e-7,1e-6,2), X=c(0,0,0))
# Run old model (uses jCO2a and jCO2p to deliver CO2 to photosynthesis SU)
run <- run_coral(time=time, env=env,
pars=replace(defpars,
c("kNPQ", "kROS", "k", "b", "sigmaCS", "sigmaCH", "jCO2a", "jHGm"),
c(60, 60, 1, 5, 0.9, 0.9, 0.32, 1)),
NPQ="fixed2")
plot_run(run)
New model
# Run new model (sets jCO2a and jCO2p to zero; CO2 is delivered by jCw)
run <- run_coral_wC(time=time, env=env,
pars=replace(defpars,
c("kNPQ", "kROS", "k", "b", "sigmaCS", "sigmaCH", "jCO2a", "jCO2p", "jHGm"),
c(60, 60, 1, 5, 0.9, 0.9, 0, 0, 1)),
NPQ="fixed2")
plot_run_wC(run)
Varying light
Old model
# Set up simulation
time <- seq(1,1230,1)
env <- init_env(time=time, L=c(10,25,2), N=c(1e-7,1e-7,2), X=c(0,0,0))
# Run old model (uses jCO2a and jCO2p to deliver CO2 to photosynthesis SU)
run <- run_coral(time=time, env=env,
pars=replace(defpars,
c("kNPQ", "kROS", "k", "b", "sigmaCS", "sigmaCH", "jCO2a"),
c(60, 60, 1, 5, 0.9, 0.9, 0.32)),
NPQ="fixed2")
plot_run(run)
New model
# Run new model (sets jCO2a and jCO2p to zero; CO2 is delivered by jCw)
run <- run_coral_wC(time=time, env=env,
pars=replace(defpars,
c("kNPQ", "kROS", "k", "b", "sigmaCS", "sigmaCH", "jCO2a", "jCO2p"),
c(60, 60, 1, 5, 0.9, 0.9, 0, 0)),
NPQ="fixed2")
plot_run_wC(run)
Light stress
old model: Increasing Light in a non-feeding host
# Set up simulation
env <- init_env(time=time, L=c(2,60,0), N=c(1e-7,1e-7,0), X=c(0,0,0))
# Run new model (sets jCO2a and jCO2p to zero; CO2 is delivered by jCw)
run <- run_coral(time=time, env=env,
pars=replace(defpars,
c("kNPQ", "kROS", "k", "b", "sigmaCS", "sigmaCH", "jCO2a"),
c(60, 60, 1, 5, 0.9, 0.9, 0.32)),
NPQ="fixed2")
plot_run(run)
new model: Increasing Light in a non-feeding host
# Set up simulation
env <- init_env(time=time, L=c(2,60,0), N=c(1e-7,1e-7,0), X=c(0,0,0))
# Run new model (sets jCO2a and jCO2p to zero; CO2 is delivered by jCw)
run <- run_coral_wC(time=time, env=env,
pars=replace(defpars,
c("kNPQ", "kROS", "k", "b", "sigmaCS", "sigmaCH", "jCO2a", "jCO2p"),
c(60, 60, 1, 5, 0.9, 0.9, 0, 0)),
NPQ="fixed2")
plot_run_wC(run)
new model: Increasing Light in a feeding host
# Set up simulation
env <- init_env(time=time, L=c(2,60,0), N=c(1e-7,1e-7,0), X=c(5e-6,5e-6,0))
# Run new model (sets jCO2a and jCO2p to zero; CO2 is delivered by jCw)
run <- run_coral_wC(time=time, env=env,
pars=replace(defpars,
c("kNPQ", "kROS", "k", "b", "sigmaCS", "sigmaCH", "jCO2a", "jCO2p"),
c(60, 60, 1, 5, 0.9, 0.9, 0, 0)),
NPQ="fixed2")
plot_run_wC(run)
Light stress with old model
# Set up simulation
time <- seq(1,1230,1)
env <- init_env(time=time, L=c(10,40,2), N=c(1e-7,1e-7,2), X=c(0,0,0))
# Run old model (uses jCO2a and jCO2p to deliver CO2 to photosynthesis SU)
run <- run_coral(time=time, env=env,
pars=replace(defpars,
c("kNPQ", "kROS", "k", "b", "sigmaCS", "sigmaCH", "jCO2a"),
c(60, 60, 1, 5, 0.9, 0.9, 0.32)),
NPQ="fixed2")
plot_run(run)
Light stress with new model
# Set up simulation
time <- seq(1,1230,1)
env <- init_env(time=time, L=c(10,40,2), N=c(1e-7,1e-7,2), X=c(0,0,0))
# Run new model (sets jCO2a and jCO2p to zero; CO2 is delivered by jCw)
run <- run_coral_wC(time=time, env=env, kCO2=2,
pars=replace(defpars,
c("kNPQ", "kROS", "k", "b", "sigmaCS", "sigmaCH", "jCO2a", "jCO2p"),
c(60, 60, 1, 5, 0.9, 0.9, 0, 0)),
NPQ="fixed2")
plot_run_wC(run)
Examine sources of CO2 for photosynthesis in new model
# Run new model (sets jCO2a and jCO2p to zero; CO2 is delivered by jCw and metabolic CO2)
run <- run_coral_wC(time=time, env=env,
pars=replace(defpars,
c("kNPQ", "kROS", "k", "b", "sigmaCS", "sigmaCH", "jCO2a", "jCO2p"),
c(60, 60, 1, 5, 0.9, 0.9, 0, 0)),
NPQ="fixed2")
par(mfrow=c(4,1), mar=c(2,2,1,1))
# Plot light
plot_L(run)
# Plot symbiont to host ratio
plot_sh(run)
# Plot symbiont-specific CO2 delivery sources
with(run, {
totalCO2 <- (H$jCw+H$rCH)*H$H/S$S + S$rCS
plot(time, totalCO2, type="l", col="gray")
lines(time, H$jCw*H$H/S$S, lty=1, col="black")
lines(time, H$rCH*H$H/S$S, lty=2, col="brown")
lines(time, S$rCS, lty=2, col="green")
})
# Plot proportions of each CO2 source
with(run, {
totalCO2 <- (H$jCw+H$rCH)*H$H/S$S + S$rCS
plot(time, H$jCw*H$H/S$S/totalCO2, type="l", ylim=c(0,1))
lines(time, H$rCH*H$H/S$S/totalCO2, lty=2, col="brown")
lines(time, S$rCS/totalCO2, lty=2, col="green")
})
Light stress with new model, +feeding
# Set up simulation
env <- init_env(time=time, L=c(10,40,2), N=c(1e-7,1e-7,2), X=c(1e-6,1e-6,0))
# Run new model (sets jCO2a and jCO2p to zero; CO2 is delivered by jCw)
run <- run_coral_wC(time=time, env=env,
pars=replace(defpars,
c("kNPQ", "kROS", "k", "b", "sigmaCS", "sigmaCH", "jCO2a", "jCO2p"),
c(60, 60, 1, 5, 0.9, 0.9, 0, 0)),
NPQ="fixed2")
plot_run_wC(run)
Light stress with new model, +feeding, reduced DIN
# Set up simulation
env <- init_env(time=time, L=c(10,40,2), N=c(0.5e-7,0.5e-7,2), X=c(1e-6,1e-6,0))
# Run new model (sets jCO2a and jCO2p to zero; CO2 is delivered by jCw)
run <- run_coral_wC(time=time, env=env,
pars=replace(defpars,
c("kNPQ", "kROS", "k", "b", "sigmaCS", "sigmaCH", "jCO2a", "jCO2p"),
c(60, 60, 1, 5, 0.9, 0.9, 0, 0)),
NPQ="fixed2")
plot_run_wC(run)
Light stress with new model, +tolerant symbiont
# Set up simulation
env <- init_env(time=time, L=c(10,40,2), N=c(1e-7,1e-7,2), X=c(1e-6,1e-6,0))
# Run new model (sets jCO2a and jCO2p to zero; CO2 is delivered by jCw)
run <- run_coral_wC(time=time, env=env,
pars=replace(defpars,
c("kNPQ", "kROS", "k", "b", "sigmaCS", "sigmaCH", "jCO2a", "jCO2p"),
c(60, 120, 1, 5, 0.9, 0.9, 0, 0)),
NPQ="fixed2")
plot_run_wC(run)
Increasing DIN
Old model
# Set up simulation
env <- init_env(time=time, L=c(30,30,0), N=c(1e-8,4e-6,0), X=c(0,0,0))
# Run old model
run <- run_coral(time=time, env=env,
pars=replace(defpars,
c("kNPQ", "kROS", "k", "b", "sigmaCS", "sigmaCH", "jCO2a"),
c(60, 60, 1, 5, 0.9, 0.9, 0.32)),
NPQ="fixed2")
plot_run(run)
New model
# Set up simulation
env <- init_env(time=time, L=c(30,30,0), N=c(1e-8,4e-6,0), X=c(0,0,0))
# Run new model (sets jCO2a and jCO2p to zero; CO2 is delivered by jCw)
run <- run_coral_wC(time=time, env=env,
pars=replace(defpars,
c("kNPQ", "kROS", "k", "b", "sigmaCS", "sigmaCH", "jCO2a", "jCO2p"),
c(60, 60, 1, 5, 0.9, 0.9, 0, 0)),
NPQ="fixed2")
plot_run_wC(run)
Steady state behavior of new model
Starting biomass = 0.01
# Set time and parameters
time <- seq(1,200,0.5)
defpars <- def_pars()
defpars <- replace(defpars, "jHGm", 1)
pars <- replace(defpars, c("kNPQ", "kROS", "k", "b", "sigmaCS", "sigmaCH", "jCO2a", "jCO2p", "initS"),
c(60, 60, 1, 5, 0.9, 0.9, 0, 0, 0.01))
# Set values of L and N at which to get steady state values
at <- expand.grid(L=seq(from=0, to=60, length.out=21),
N=seq(from=0, to=4e-6, length.out=21))
# Run to steady state and plot results
ssdat <- run_steady_states(time=time, pars=pars, at=at, food=0)
plot_steady_states(ssdat)
Starting biomass = 0.1
# Set time and parameters
time <- seq(1,200,0.5)
defpars <- def_pars()
defpars <- replace(defpars, "jHGm", 1)
pars <- replace(defpars, c("kNPQ", "kROS", "k", "b", "sigmaCS", "sigmaCH", "jCO2a", "jCO2p", "initS"),
c(60, 60, 1, 5, 0.9, 0.9, 0, 0, 0.1))
# Set values of L and N at which to get steady state values
at <- expand.grid(L=seq(from=0, to=60, length.out=21),
N=seq(from=0, to=4e-6, length.out=21))
# Run to steady state and plot results
ssdat <- run_steady_states(time=time, pars=pars, at=at, food=0)
plot_steady_states(ssdat)
Starting biomass = 0.5
# Set time and parameters
time <- seq(1,200,0.5)
defpars <- def_pars()
defpars <- replace(defpars, "jHGm", 1)
pars <- replace(defpars, c("kNPQ", "kROS", "k", "b", "sigmaCS", "sigmaCH", "jCO2a", "jCO2p", "initS"),
c(60, 60, 1, 5, 0.9, 0.9, 0, 0, 0.5))
# Set values of L and N at which to get steady state values
at <- expand.grid(L=seq(from=0, to=60, length.out=21),
N=seq(from=0, to=4e-6, length.out=21))
# Run to steady state and plot results
ssdat <- run_steady_states(time=time, pars=pars, at=at, food=0)
plot_steady_states(ssdat)
Starting symbiont biomass = 1
# Set time and parameters
time <- seq(1,200,0.5)
pars <- replace(defpars, c("kNPQ", "kROS", "k", "b", "sigmaCS", "sigmaCH", "jCO2a", "jCO2p", "initS"),
c(60, 60, 1, 5, 0.9, 0.9, 0, 0, 1))
# Set values of L and N at which to get steady state values
at <- expand.grid(L=seq(from=0, to=60, length.out=21),
N=seq(from=0, to=4e-6, length.out=21))
# Run to steady state and plot results
ssdat <- run_steady_states(time=time, pars=pars, at=at, food=0)
plot_steady_states(ssdat)
+ feeding
# Set time and parameters
time <- seq(1,100,1)
pars <- replace(defpars, c("kNPQ", "kROS", "k", "b", "sigmaCS", "sigmaCH", "jCO2a", "jCO2p"),
c(60, 60, 1, 5, 0.9, 0.9, 0, 0))
# Set values of L and N at which to get steady state values
at <- expand.grid(L=seq(from=0, to=60, length.out=41),
N=seq(from=0, to=4e-6, length.out=41))
# Run to steady state and plot results
ssdat <- run_steady_states(time=time, pars=pars, at=at, food=5e-6)
plot_steady_states(ssdat)
+ feeding + tolerant symbiont
# Set time and parameters
time <- seq(1,100,1)
pars <- replace(defpars, c("kNPQ", "kROS", "k", "b", "sigmaCS", "sigmaCH", "jCO2a", "jCO2p"),
c(60, 200, 1, 5, 0.9, 0.9, 0, 0))
# Set values of L and N at which to get steady state values
at <- expand.grid(L=seq(from=0, to=60, length.out=41),
N=seq(from=0, to=4e-6, length.out=41))
# Run to steady state and plot results
ssdat <- run_steady_states(time=time, pars=pars, at=at, food=5e-6)
plot_steady_states(ssdat)
Stability at high nitrogen…?
When system starts at 2 µM DIN and increases to 4 µM DIN, switch is btw. 0.15 and 0.2 initS.
# Set up simulation
env <- init_env(time=time, L=c(10,10,0), N=c(1e-6,4e-6,0), X=c(0,0,0))
# Run new model (sets jCO2a and jCO2p to zero; CO2 is delivered by jCw)
run <- run_coral_wC(time=time, env=env,
pars=replace(defpars,
c("kNPQ", "kROS", "k", "b", "sigmaCS", "sigmaCH", "jCO2a", "jCO2p", "initS"),
c(60, 60, 1, 5, 0.9, 0.9, 0, 0, 0.15)),
NPQ="fixed2")
plot_run_wC(run)
# Set up simulation
env <- init_env(time=time, L=c(10,10,0), N=c(1e-6,4e-6,0), X=c(0,0,0))
# Run new model (sets jCO2a and jCO2p to zero; CO2 is delivered by jCw)
run <- run_coral_wC(time=time, env=env,
pars=replace(defpars,
c("kNPQ", "kROS", "k", "b", "sigmaCS", "sigmaCH", "jCO2a", "jCO2p", "initS"),
c(60, 60, 1, 5, 0.9, 0.9, 0, 0, 0.2)),
NPQ="fixed2")
plot_run_wC(run)
When system starts at 2 µM DIN and increases to 4 µM DIN, switch is btw. 0.2 and 0.3 initS.
# Set up simulation
env <- init_env(time=time, L=c(10,10,0), N=c(2e-6,4e-6,0), X=c(0,0,0))
# Run new model (sets jCO2a and jCO2p to zero; CO2 is delivered by jCw)
run <- run_coral_wC(time=time, env=env, kCO2=5,
pars=replace(defpars,
c("kNPQ", "kROS", "k", "b", "sigmaCS", "sigmaCH", "jCO2a", "jCO2p", "jHGm", "initS"),
c(60, 60, 1, 5, 0.9, 0.9, 0, 0, 1, 0.2)),
NPQ="fixed2")
plot_run_wC(run)
# Set up simulation
env <- init_env(time=time, L=c(9,9,0), N=c(2e-6,4e-6,0), X=c(0,0,0))
# Run new model (sets jCO2a and jCO2p to zero; CO2 is delivered by jCw)
run <- run_coral_wC(time=time, env=env, kCO2=5,
pars=replace(defpars,
c("kNPQ", "kROS", "k", "b", "sigmaCS", "sigmaCH", "jCO2a", "jCO2p", "jHGm", "initS"),
c(60, 60, 1, 5, 0.9, 0.9, 0, 0, 1, 0.3)),
NPQ="fixed2")
plot_run_wC(run)
# Set up simulation
env <- init_env(time=time, L=c(10,10,0), N=c(4e-6,4e-6,0), X=c(0,0,0))
# Run new model (sets jCO2a and jCO2p to zero; CO2 is delivered by jCw)
run <- run_coral_wC(time=time, env=env, kCO2=5,
pars=replace(defpars,
c("kNPQ", "kROS", "k", "b", "sigmaCS", "sigmaCH", "jCO2a", "jCO2p", "jHGm", "initS"),
c(60, 60, 1, 5, 0.9, 0.9, 0, 0, 1, 0.5)),
NPQ="fixed2")
plot_run_wC(run)
# Set up simulation
env <- init_env(time=time, L=c(10,10,0), N=c(4e-6,4e-6,0), X=c(0,0,0))
# Run new model (sets jCO2a and jCO2p to zero; CO2 is delivered by jCw)
run <- run_coral_wC(time=time, env=env, kCO2=5,
pars=replace(defpars,
c("kNPQ", "kROS", "k", "b", "sigmaCS", "sigmaCH", "jCO2a", "jCO2p", "jHGm", "initS"),
c(60, 60, 1, 5, 0.9, 0.9, 0, 0, 1, 0.01)),
NPQ="fixed2")
plot_run_wC(run)
with(run, {
plot(S$S/H$H, S$jSG)
})